home *** CD-ROM | disk | FTP | other *** search
/ Amiga Format CD 24 / Amiga Format AFCD24 (Feb 1998, Issue 108).iso / -in_the_mag- / emulation / adfs / coma.adf / BinaryGCD.asm < prev    next >
Assembly Source File  |  1998-01-04  |  6KB  |  209 lines
  1. * 68K implementation of Donald Knuth's binary Greatest Common Divisor algorithm
  2. * The Art of Computer Programming, Volume 2, Seminumerical Algorithms, page 321
  3. * Discovered by J Stein in 1961, published in J Comp Phys 1, page 397-405, 1967
  4. * Written in generic 68000 code and tested on 68060 by Simon N Goodwin 31/12/97
  5. *
  6. * To do: register equates, tests on more CPUs
  7. *
  8. * Register parameters
  9. *
  10. * Input:   D0 and D1, POSITIVE non-zero words to be tested (U and V)
  11. * Result:  D0, greatest common divisor of the inputs (R)
  12. *
  13. * Variables
  14. *
  15. *    K          D2       ; Shift count for common powers of 2
  16. *     T         D3       ; Contender to replace larger of U and V
  17. *     S    D4    ; Scratch
  18. *
  19. * All values are treated as 32 bit signed long words
  20. *
  21. *******************************************************************************
  22. *
  23. * Test bed - returns number of pairs with common factor >1, should be 40% ish
  24. * If any 'optimisation' changes this number, you must have broken something!
  25. *
  26. GEN2    EQU    1    ; 1 for 2 random calls per GCF, 0 for one
  27. LOOPS    EQU    10000
  28. *
  29. Test    move.l    #LOOPS,d7    ; Loop count
  30.     moveq    #0,d5    ; Factor count
  31.     bra.s    loop    ; Just in case CNOP doesn't
  32.     nop
  33. *
  34. * Generate two random numbers. The generator does produce two at once
  35. * in D0 and D1 so only one call is really needed - but D0 is not very
  36. * random, hence 2 distinct calls for each test are preferable (GEN2).
  37. *
  38.     cnop    0,8    ; Align loop start
  39.     
  40. loop    jsr    LongRand60
  41.     move.l    d1,d6
  42.     lsr.l    #1,d6    ; Keep it positive
  43.     beq.s    loop    ; Avoid zero
  44.  
  45. * Set GEN2 to 0 to use high and low words of one 64 bit random number
  46. * instead of two 32 bit values generated separately (with GEN2 EQU 1)
  47.  
  48.     ifne    GEN2    
  49.  
  50.     jsr    LongRand60  ; Get another rather random set
  51.     lsr.l    #1,d1    ; Stay positive, most random word
  52.     beq.s    loop    ; Unlikely but fatal!
  53.     move.l    d6,d0    
  54.  
  55.     else
  56.  
  57.     lsr.l    #1,d0    ; Stay positive
  58.     beq.s    loop    ; Unlikely but fatal!
  59.     move.l    d6,d1
  60.  
  61.     endc
  62.     
  63. Try1    jsr    GCD_68000    ; Do the real work
  64.  
  65.     subq.l    #1,d0
  66.     beq.s    no_divisor
  67.  
  68.     addq.l    #1,d5    ; Count one more success
  69.  
  70. no_divisor    subq.l    #1,d7    ; Keep trying
  71.     bgt.s    loop    
  72.  
  73. Done    
  74. *    moveq    #0,d0
  75.     move.l    d5,d0    ; Return success count
  76.     rts    
  77. *
  78. * Long word polynomial pseudo-random number generator
  79. * This is optimised for 68060s but runs on others too
  80. *
  81.     cnop    0,8    ; Align subroutine start
  82. *    
  83. LongRand60    move.l    Rand60+4,d0 ; Soep
  84.     move.l    Rand60,d1    ; Poep
  85.     andi.b    #14,d0    ; Soep
  86.     ori.b    #32,d0    ; Poep
  87.     move.l    d1,d3    ; Soep
  88.     move.l     d0,d2    ; Poep
  89.     add.l      d2,d2    ; Soep (bypass in)
  90.     addx.l     d3,d3    ; Poep only
  91.     add.l      d2,d0    ; Poep
  92.     moveq      #16,d4    ; Soep - rotation factor
  93.     addx.l     d3,d1    ; Poep only
  94.     ror.l      d4,d3    ; Poep
  95.     ror.l      d4,d2    ; Soep
  96.     move.w     d2,d3    ; Poep
  97.     clr.w      d2    ; Soep
  98.     add.l      d2,d0    ; Poep
  99.     move.l     d0,Rand60+4  ; Soep (bypass out)
  100.     addx.l     d3,d1    ; Poep only
  101.     move.l     d1,Rand60    ; Poep
  102.     rts
  103. *    
  104.     cnop    0,8    ; Long align data (or better)
  105. Rand60    dc.l    0,0    ; Seed
  106. *
  107. * 68060 time for a million tests with generic GCD code: 10.68 seconds
  108. * RTS at start of GCD_68000 = without generic GCD code: 1.86 seconds
  109. *
  110. * Thus the total loop overhead including random number generation is
  111. * less than 20 per cent - clearly GCD computation time predominates.
  112. * All times are approximate as I'm testing an active (but not busy)
  113. * system without resorting to the use of FORBID/DISABLE pairs, etc.
  114. *
  115. * The number of calls to the random function may be halved by using
  116. * both halves of the 64 bit random number as two 31 bit numbers:
  117. *
  118. * 68060 time for a million tests with generic GCD code: 10.05 seconds
  119. * RTS at start of GCD_68000 = without generic GCD code: 1.15 seconds
  120. *
  121. * Thus the test rig is contributing getting on for as much time as
  122. * one pass through the randomising function (which seems to be taking
  123. * around two thirds of a second or perhaps 30 T states, suggesting that
  124. * the subroutine call overhead may be the next thing usefully trimmed.
  125. * Just adding some NOPs saved three to four per cent of total run time!
  126. *
  127. * 68060 GEN2 set, 1 million tests with generic GCD code: 10.28 seconds
  128. * 68060 GEN2 unset, million tests with generic GCD code:  9.74 seconds
  129. *
  130. *******************************************************************************
  131. *
  132. * Generic 68000 version - first and simplest version
  133. *
  134. * This is an attempt at a straight-forward solution of the problem,
  135. * but a few 68000-specific performance peculiarities are indulged.
  136. *
  137. * Three days trying to produce a faster version for the 68060 failed,
  138. * including BTST #0s and a restoring version of the B3 loop, so this
  139. * is the definitive 68K version, unless I think of something better!
  140. *
  141. * IF ANY INPUT IS NEGATIVE OR ZERO THIS ALGORITHM MAY NOT TERMINATE!
  142. *
  143.     cnop    0,8    ; Align branch target
  144. *    
  145. GCD_68000
  146. *
  147. * First stage, count and discard low-order zero bits in both operands
  148. *
  149. B1    moveq    #0,d2    ; Initial shift count
  150. S1    move.l    d0,d4    ; Combine operands U an V
  151.     or.l    d1,d4
  152.     and.w    #1,d4    ; Test low order bit 
  153.             ; (.b is no faster and .l is longer)
  154.     bne.s    B2    ; Exit unless both U and V are even
  155.     
  156.     asr.l    #1,d0    ; Shift out zero bit
  157.     asr.l    #1,d1    ; Do the same for both operands
  158.     addq.l    #1,d2    ; Count one power of two in factor
  159.     bra.s    S1    ; Try again
  160. *
  161. * Second stage, initialise T to (U & 1) ? -V : U
  162. *
  163. B2    moveq    #1,d4    ; Mask for least significant bit
  164.     and.w    d0,d4    ; Test bit 0 of U (BTST #0 may be slower)
  165.     beq.s    U_even    ; Continue only if U is odd
  166.     
  167.     move.l    d1,d3    ; Set T to V
  168.     neg.l    d3    ; Set T to -V
  169.     bra.s    B4    ; Some CPUs favour TRAPF.L ($51FB)
  170.  
  171. U_even    move.l    d0,d3    ; Set T to U
  172. *
  173. * Third stage, T is now non-zero and even. Halve it.
  174. *
  175. B3    asr.l    #1,d3
  176. *
  177. * Fourth stage, See if T is still even
  178. *
  179. B4    moveq    #1,d4    ; Bit 0 mask
  180.     and.w    d3,d4    ; Test low order bit (word, see above)
  181.     beq.s    B3    ; Iterate until T becomes odd
  182. *
  183. * Fifth stage, Reset MAX(U,V) so larger of U and V is replaced by ABS(T)
  184. * T > 0 ? U := T : V := -T
  185. *
  186. B5    move.l    d3,d4    ; Copy T to scratch, implicitly testing it
  187.     bgt.s    Replace_U
  188.     
  189.     neg.l    d4    ; Generate -T
  190.     move.l    d4,d1    ; Set V to -T
  191.     move.l    d0,d3    ; Set T to U in preparation for stage 6
  192.     
  193. Replace_U    move.l    d3,d0    ; Set U to T - not worth skipping!
  194. *
  195. * Sixth stage, Set T to U-V and iterate till T=0
  196. *
  197. B6    sub.l    d1,d3    ; Set T to U-V
  198.     bne.s    B3    ; Iterate while T<>0
  199. *
  200. * Epilogue, result is U * 2 ^ K
  201. *
  202.     lsl.l    d2,d0    ; Scale the result in D0
  203.     rts
  204. *
  205. *******************************************************************************
  206. *
  207. *    dc.w    $51FA    ; Test TRAPF.W, skip word; $51FB skips long
  208.  
  209.     end